#!/usr/bin/perl -w
use strict;
use DBI;
use DBD::Pg;
use FindBin;
use lib ("$FindBin::Bin/..", "/net/cpp-group/Leo/bin");
use uniquify;
use run_cmd;
use die_error;

use db_parameters;
use ortho_species;

print <<"HEADLINE";
888888888888888888888888888888888888888888888888888888888888888888888

    table_closest_same_species

		Get the same_species protein or gene with the lowest dS

888888888888888888888888888888888888888888888888888888888888888888888
HEADLINE




#
#	Connect to genome database
#
my $dbh = connect_to_panda();
use get_psql_time;
my $start_time = get_psql_time($dbh);
#my @db_param = get_db_parameters();


	#
	#	make map from prot_ids and gene_ids
	#
	{
		open_or_die (*MAP, "$dir_orthologs_output/prot_to_gene_ids.map");
		open_or_die (*NEWMAP, ">$dir_pipeline_temp/prot_to_prot_gene_ids.map");
		while (<MAP>)
		{
			my ($prot_id,  $gene_id,  $species) = split /\t/;
			print NEWMAP $prot_id. "\t". $prot_id.'<TAB>'.$gene_id. "\n";
		}
	}

	#
	#	map ds to prot_id???gene_id
	#
	print STDERR "\tMap ds from prot_id to gene_id\n";
	my $cmd = <<"    CMD";
	seq_pairs_substitute_identifiers
		--input_file $dir_orthologs_output/ds_res.both_filtered
		--map_identifiers $dir_pipeline_temp/prot_to_prot_gene_ids.map
		--verbose
		>
		$dir_pipeline_temp/prot_gene_ids.ds
    CMD
	run_cmd($cmd);


unlink("$dir_pipeline_temp/same_species.ds");
for my $species($ortho_gene1,$ortho_gene2)
{
	
	#
	# Get lowest dS
	#
	print STDERR "\n\tGet dS for all pairs from the same species...\n";
	my $cmd = <<"    CMD";
	seq_pairs_filter_by_regex
							--regex "($species)"
                            --patterns "$species"
							--verbose 
							--verbose 
                            --input_file $dir_pipeline_temp/prot_gene_ids.ds
							>>
							$dir_pipeline_temp/same_species.ds
    CMD
	run_cmd($cmd);
}


	print STDERR "\n\tFor each protein, get pairs with the lowest dS...\n";
	$cmd = <<"    CMD";
	seq_pairs_lowest
							--regex "<TAB>(.*)"
							--reciprocals
							--verbose
                            --input_file $dir_pipeline_temp/same_species.ds
							 > $dir_pipeline_temp/closest_same_species.ds
    CMD
	run_cmd($cmd);


	open DS, "$dir_pipeline_temp/closest_same_species.ds" 
							or die "Error:\n\tCould not open $dir_pipeline_temp/closest_same_species.ds\n";
	print STDERR "$dir_pipeline_temp/closest_same_species.ds\n";
							
	print STDERR "\n\tEntering into panda...\n";
	$dbh->do("CREATE TEMP TABLE t1 (prot_id TEXT, gene_id TEXT, same_species_prot_id TEXT, same_species_gene_id TEXT, ds REAL);");
	$dbh->do("COPY t1".
					"(prot_id, gene_id, same_species_prot_id, same_species_gene_id, dS) FROM STDIN") or die $dbh->errstr.".";;
	while (<DS>)
	{
		s/<TAB>/\t/g or die_error("Ill formed line \n[" . $_ . "]");
		$dbh->pg_putline($_);
	}
	$dbh->pg_endcopy();

	print STDERR "\n\tSetting orthologs.closest_same_species for gene_ids...\n";
	my $sql_cmd =<<"PL/SQLCMD";
	TRUNCATE TABLE orthologs.closest_same_species;
    INSERT INTO 
        orthologs.closest_same_species
		(
			gene_id,
			prot_id,
			species, 
			same_species_gene_id,
			same_species_prot_id,
			dS,
            protocol_id
		)
		SELECT
    			DISTINCT ON (gene_id)
    			gene_id,
				prot_id,
    			species, 
    			same_species_gene_id,
    			same_species_prot_id,
    			dS,
    			$curr_protocol_id
    		FROM
    			t1 NATURAL JOIN 
    			ens_id;

	VACUUM ANALYSE orthologs.closest_same_species;
PL/SQLCMD
	$dbh->do($sql_cmd);
	


	
	print STDERR "\tDone\n";
	print STDERR "\tCompleted\n";

	unlink("$dir_pipeline_temp/prot_gene_ids.ds");
	unlink("$dir_pipeline_temp/prot_to_prot_gene_ids.map");
	unlink("$dir_pipeline_temp/same_species.ds");
	unlink("$dir_pipeline_temp/closest_same_species.ds");
	log_pipeline_stage($dbh, $curr_protocol_id, 32, 'table closest_same_species', $start_time);
	$dbh->disconnect();

